{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.fftpack import fft,ifft\n",
    "import pandas as pd\n",
    "from scipy.optimize import minimize\n",
    "from iminuit import Minuit\n",
    "from scipy.optimize import basinhopping\n",
    "from scipy.optimize import curve_fit\n",
    "import emcee\n",
    "from pprint import pprint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 似然函数 p; D = -2 ln p\n",
    "\n",
    "def twi_minus_loglikelihood(A,f_b,alpha_H,poisson):\n",
    "    alpha_L = 1.0\n",
    "    \n",
    "    perdata00 = pd.read_csv(\"perlist00.csv\")\n",
    "    f = perdata00['f']\n",
    "    per = perdata00['per']\n",
    "            \n",
    "    model = []\n",
    "    f_length = len(f)\n",
    "    for i in range(f_length):\n",
    "        model.append(((f[i]**(-alpha_L))/(1+(f[i]/f_b)**(alpha_H-alpha_L)))*A+poisson)\n",
    "     \n",
    "    \n",
    "    length = len(perdata00)\n",
    "    minus_log_p = 0\n",
    "    for i in range(length):\n",
    "        minus_log_p += (per[i]/model[i]+math.log(model[i]))\n",
    "    \n",
    "    \n",
    "    D = 2*minus_log_p\n",
    "    print (D)\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "m=Minuit(twi_minus_loglikelihood,A=0.005,f_b=1.7E-4,alpha_H=3.8,poisson=0.8,\n",
    "         error_A=0.0001,error_f_b=1.0E-5,error_alpha_H=0.01,error_poisson=0.01,\n",
    "         limit_A=(0.001,0.01), limit_f_b=(1.0E-4,1.0E-3),limit_alpha_H=(2.0,5.0),limit_poisson=(0,1),\n",
    "         errordef=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "326.25985765080605\n",
      "326.25696216125806\n",
      "326.26275464750574\n",
      "326.2309704812444\n",
      "326.2888955356907\n",
      "326.238286435313\n",
      "326.28146204725385\n",
      "326.2131195892432\n",
      "326.3067517618009\n",
      "326.2609006688678\n",
      "326.2588147016205\n",
      "326.27029092772545\n",
      "326.24943126153454\n",
      "326.2917124243864\n",
      "326.22799805916424\n",
      "326.44108150244955\n",
      "326.07847774976824\n",
      "326.07847774976824\n",
      "325.3514131255773\n",
      "323.5232116486788\n",
      "317.9555109343111\n",
      "300.7263168898887\n",
      "253.79464117376818\n",
      "751.095161926842\n",
      "277.0340613945565\n",
      "247.72698665139416\n",
      "245.7874266494233\n",
      "245.50019523996477\n",
      "245.4498046393886\n",
      "245.55077347742431\n",
      "245.41930370624232\n",
      "245.58128110640308\n",
      "245.5989888345741\n",
      "245.4015114139126\n",
      "245.48318621340297\n",
      "245.518366521132\n",
      "245.49466041266905\n",
      "245.50584760750957\n",
      "249.69331675391456\n",
      "241.37763133717965\n",
      "233.9799399166395\n",
      "226.9852859623153\n",
      "226.98474706986406\n",
      "226.98588421868226\n",
      "226.98413677165286\n",
      "226.986485431202\n",
      "226.98361352314745\n",
      "226.98706711330286\n",
      "226.98805253564376\n",
      "226.98292544823917\n",
      "226.98666683067827\n",
      "226.98401380593293\n",
      "226.9985711625779\n",
      "226.9720893773363\n",
      "225.9299535570379\n",
      "225.87925169129122\n",
      "225.88160426766166\n",
      "225.87700588253023\n",
      "225.88095782721643\n",
      "225.8776451027204\n",
      "225.87911967073092\n",
      "225.87949782880133\n",
      "225.87888463359988\n",
      "225.8797536448299\n",
      "225.8564419523859\n",
      "225.84315322603175\n",
      "225.84349607266245\n",
      "225.84293371527502\n",
      "225.84347796018534\n",
      "225.84293522065693\n",
      "225.84336851990895\n",
      "225.84304679565298\n",
      "225.84428515344362\n",
      "225.84212712014724\n",
      "225.83328567164577\n",
      "225.83254678146318\n",
      "225.83248434225123\n",
      "225.83272292818964\n",
      "225.83278912091538\n",
      "225.83241068560775\n",
      "225.83262219265995\n",
      "225.8325804132214\n",
      "225.83259834425675\n",
      "225.8326051047113\n",
      "225.83233303464237\n",
      "225.83197433822997\n",
      "225.83223529465567\n",
      "225.83181991352546\n",
      "225.83254987952589\n",
      "225.8315066565568\n",
      "225.83173435851108\n",
      "225.83232346988868\n",
      "225.83212291253244\n",
      "225.8319338002575\n",
      "225.83177297543722\n",
      "225.83099178155337\n",
      "225.82920529267997\n",
      "225.82519430186065\n",
      "225.82259273272928\n",
      "225.82251303043088\n",
      "225.82278145602598\n",
      "225.82233697381545\n",
      "225.82267835263326\n",
      "225.82245565384991\n",
      "225.82284792514304\n",
      "225.82230024205705\n",
      "225.8221599786033\n",
      "225.8229749691104\n",
      "225.8189679415461\n",
      "225.81889649689938\n",
      "225.81912311570608\n",
      "225.81877488700425\n",
      "225.81913303355392\n",
      "225.818765236944\n",
      "225.81878402816636\n",
      "225.81912274599705\n",
      "225.81899641542606\n",
      "225.81890404838538\n",
      "225.81874023099374\n",
      "225.8188107939156\n",
      "225.81877802171732\n",
      "225.81880610668944\n",
      "225.81878252052059\n",
      "225.81879037800064\n",
      "225.81879841884665\n",
      "225.81881454958426\n",
      "225.81877409446082\n",
      "225.81874023099374\n",
      "225.8188107939156\n",
      "225.81877802171732\n",
      "225.81880610668944\n",
      "225.81878252052059\n",
      "225.81879037800064\n",
      "225.81879841884665\n",
      "225.81881454958426\n",
      "225.81877409446082\n",
      "225.81874569458833\n",
      "225.81873910153763\n",
      "225.81874475246207\n",
      "225.8187400361387\n",
      "225.81874159430564\n",
      "225.8187432010799\n",
      "225.8187464454386\n",
      "225.8187383438296\n",
      "225.818978241607\n",
      "225.81878026043063\n",
      "225.81890511888713\n",
      "225.81875811933256\n",
      "225.81890564130026\n",
      "225.81882373683584\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td title=\"Minimum value of function\">FCN = 225.81874023099374</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 150</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 5.131682437190387e-06</td>\n",
       "        <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 1e-05</td>\n",
       "        <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
       "        UP = 1.0</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<table>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
       "        <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
       "        <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
       "        <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
       "        <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
       "        <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
       "        <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#iwzgbEztKP').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>A</td>\n",
       "        <td>0.00593157</td>\n",
       "        <td>0.00676127</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.001</td>\n",
       "        <td>0.01</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>f_b</td>\n",
       "        <td>0.00016548</td>\n",
       "        <td>0.000183371</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.0001</td>\n",
       "        <td>0.001</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>alpha_H</td>\n",
       "        <td>2.58801</td>\n",
       "        <td>0.750826</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>2</td>\n",
       "        <td>5</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>3</td>\n",
       "        <td>poisson</td>\n",
       "        <td>0.408119</td>\n",
       "        <td>0.0296436</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"iwzgbEztKP\" style=\"display:none;\">\n",
       "<textarea rows=\"14\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & A & 0.00593157 & 0.00676127 &  &  & 0.001 & 0.01 & No\\\\\n",
       "\\hline\n",
       "1 & $f_{b}$ & 0.00016548 & 0.000183371 &  &  & 0.0001 & 0.001 & No\\\\\n",
       "\\hline\n",
       "2 & $\\alpha_{H}$ & 2.58801 & 0.750826 &  &  & 2.0 & 5 & No\\\\\n",
       "\\hline\n",
       "3 & poisson & 0.408119 & 0.0296436 &  &  & 0.0 & 1 & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "({'fval': 225.81874023099374,\n",
       "  'edm': 5.131682437190387e-06,\n",
       "  'nfcn': 150,\n",
       "  'up': 1.0,\n",
       "  'is_valid': True,\n",
       "  'has_valid_parameters': True,\n",
       "  'has_accurate_covar': True,\n",
       "  'has_posdef_covar': True,\n",
       "  'has_made_posdef_covar': False,\n",
       "  'hesse_failed': False,\n",
       "  'has_covariance': True,\n",
       "  'is_above_max_edm': False,\n",
       "  'has_reached_call_limit': False},\n",
       " [{'number': 0,\n",
       "   'name': 'A',\n",
       "   'value': 0.005931568242742783,\n",
       "   'error': 0.006761267001087538,\n",
       "   'is_const': False,\n",
       "   'is_fixed': False,\n",
       "   'has_limits': True,\n",
       "   'has_lower_limit': True,\n",
       "   'has_upper_limit': True,\n",
       "   'lower_limit': 0.001,\n",
       "   'upper_limit': 0.01},\n",
       "  {'number': 1,\n",
       "   'name': 'f_b',\n",
       "   'value': 0.0001654802221753269,\n",
       "   'error': 0.0001833711986565022,\n",
       "   'is_const': False,\n",
       "   'is_fixed': False,\n",
       "   'has_limits': True,\n",
       "   'has_lower_limit': True,\n",
       "   'has_upper_limit': True,\n",
       "   'lower_limit': 0.0001,\n",
       "   'upper_limit': 0.001},\n",
       "  {'number': 2,\n",
       "   'name': 'alpha_H',\n",
       "   'value': 2.5880147768245507,\n",
       "   'error': 0.7508261981809388,\n",
       "   'is_const': False,\n",
       "   'is_fixed': False,\n",
       "   'has_limits': True,\n",
       "   'has_lower_limit': True,\n",
       "   'has_upper_limit': True,\n",
       "   'lower_limit': 2.0,\n",
       "   'upper_limit': 5.0},\n",
       "  {'number': 3,\n",
       "   'name': 'poisson',\n",
       "   'value': 0.4081187847946647,\n",
       "   'error': 0.029643609530139947,\n",
       "   'is_const': False,\n",
       "   'is_fixed': False,\n",
       "   'has_limits': True,\n",
       "   'has_lower_limit': True,\n",
       "   'has_upper_limit': True,\n",
       "   'lower_limit': 0.0,\n",
       "   'upper_limit': 1.0}])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "225.81874023099374\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#HyJbrsIqHg').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>A</td>\n",
       "        <td>0.00593157</td>\n",
       "        <td>0.00676127</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.001</td>\n",
       "        <td>0.01</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>f_b</td>\n",
       "        <td>0.00016548</td>\n",
       "        <td>0.000183371</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.0001</td>\n",
       "        <td>0.001</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>alpha_H</td>\n",
       "        <td>2.58801</td>\n",
       "        <td>0.750826</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>2</td>\n",
       "        <td>5</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>3</td>\n",
       "        <td>poisson</td>\n",
       "        <td>0.408119</td>\n",
       "        <td>0.0296436</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"HyJbrsIqHg\" style=\"display:none;\">\n",
       "<textarea rows=\"14\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & A & 0.00593157 & 0.00676127 &  &  & 0.001 & 0.01 & No\\\\\n",
       "\\hline\n",
       "1 & $f_{b}$ & 0.00016548 & 0.000183371 &  &  & 0.0001 & 0.001 & No\\\\\n",
       "\\hline\n",
       "2 & $\\alpha_{H}$ & 2.58801 & 0.750826 &  &  & 2.0 & 5 & No\\\\\n",
       "\\hline\n",
       "3 & poisson & 0.408119 & 0.0296436 &  &  & 0.0 & 1 & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pprint(m.fval)\n",
    "m.print_param()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAHqCAYAAACEBnMbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XeYXHXZ//H3nSVlIRBapGOCoYiAlA1Vio+AICQU6YiUQMRHFKyhE7FgVKqoEKQpTeAnQpAqSpMamnQMRRNECG0NPAkJyff3x9mFZUnIzmZmzpwz79d1zbU7Z2bO3FtO9pNvjZQSkiRJKo4+eRcgSZKkyhjgJEmSCsYAJ0mSVDAGOEmSpIIxwEmSJBWMAU6SJKlgDHCSJEkFY4CTJEkqGAOcJElSwSyUdwG1tvTSS6chQ4bkXYYkSdJ8PfDAA6+mlAbP73mlD3BDhgxh4sSJeZchSZI0XxHxz548r7RdqBExIiLGt7e3512KJElSVZU2wKWUJqSURg8aNCjvUiRJkqqqtAFOkiSprAxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgqmtAHOZUQkSVJZlTbAuYyIJEkqq9IGOEmSpLIywEmSJBWMAU6SJKlgDHCSJEkFY4CTJEkqGAOcJElSwZQ2wLkOnCRJKqvSBjjXgZMkSWVV2gAnSZJUVga4Khg7Nu8KJElSMzHASZIkFYwBTpIkqWAMcJIkSQVjgJMkSSoYA5wkSVLBGOAkSZIKxgAnSZJUMKUNcG6lJUmSyqq0Ac6ttCRJUlmVNsBJkiSVlQFOkiSpYAxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGACdJklQwBjhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBWOAkyRJKhgDnCRJUsGUNsBFxIiIGN/e3p53KZIkSVVV2gCXUpqQUho9aNCgvEuRJEmqqtIGOEmSpLIywEmSJBWMAW5BjR/P7lfsDo8+mnclkiSpSRjgFtS++/LiChvBNtvAF78IjzySd0WSJKnkDHALapFFuGvT78Bzz8Fmm8F228Euu8BDD+VdmSRJKikDXLUsvDB861vw7LOw5Zawww6w007wwAN5VyZJkkrGAFdtCy8MRxyRBbnPfQ5GjoQRI2DixLwrkyRJJWGAq5XWVvjGN7Ig9/nPw847Z61y992Xd2WSJKngDHC1NmAAHHZYFuR22AF22w223x7uuSfvyiRJUkEZ4Oqlf3/43/+Ff/wjGxu3555Zy9xdd+VdmSRJKhgDXL317w+HHpoFuS9+EfbZJ1uC5M47865MkiQVhAEuL/36wejR8MwzWWvcl7+cTXq4/fa8K5MkSQ3OAJe3fv3g4IPh6adh333hwAPhs5+FW2/NuzJJktSgDHCNom9fOOggeOop2H//LNRtuSX85S+QUt7VSZKkBmKAazR9+8IBB2RBbtSobLzcFlvAn/9skJMkSYABrnEttFA2Lu6JJ+ArX8mWIvnMZ+CmmwxykiQ1OQNco1toIfjSl+Dxx+FrX4PDD4dNN4UbbjDISZLUpAxwRdHSki058thjWYj79rdhk03guusMcpIkNRkDXNG0tMBee8Gjj8K3vgXf+x5stBFce61BTpKkJmGAK6o+fWCPPeDvf89C3NFHw/DhcM01BjlJkkquUAEuInaOiHMi4uqI2DbvehpCnz7Z/qoPP5yFuOOPhw02gD/+0SAnSVJJ5R7gIuK8iHglIh7rdny7iHg6IiZFxJEAKaU/ppQOAQ4A9syh3MbVpw/suis8+GAW4k48EdZbD/7wB5gzJ+/qJElSFeUe4IALgO26HoiIFuCXwPbAmsDeEbFml6cc2/G4uuvTB3beGR54IAtxP/pRFuSuvNIgJ0lSSeQe4FJKtwOvdzu8ITAppfRcSmkmcBmwU2TGAdenlB6sd62FEgEjR8LEiVmIGzcOPv1puPxyg5wkSQWXe4CbhxWAyV3uT+k49nVga2C3iDh0Xi+OiNERMTEiJk6dOrW2lTa6CNhxR7jvvizEnXwyrL02XHYZzJ6dd3WSJKkXGjXAxVyOpZTSGSmlDVJKh6aUzprXi1NK41NKbSmltsGDB9ewzAKJgC98Ae65Jwtxp5+eBblLLjHISZJUMI0a4KYAK3W5vyLw75xqKZcI2G47uOsuOO00+OUv4VOfgosugnffzbs6SZLUA40a4O4HVo2IoRHRD9gLuKaSE0TEiIgY397eXpMCCy8Ctt0W7rwTzjwTzj4b1lwTfvtbg5wkSQ0u9wAXEZcCdwOrR8SUiBiVUnoXOAy4EXgSuDyl9Hgl500pTUgpjR40aFD1iy6TCNh6a7j9djjrLDj3XPjkJ+HCCw1ykiQ1qNwDXEpp75TScimlvimlFVNK53Ycvy6ltFpK6RMppR/lXWfpRcD//A/cdhuccw5ccAGssQacfz7MmpV3dZIkqYvcA5wa0FZbwV//Cuedl42NW331rGXOICdJUkMwwGnettgCbrklGxd32WWw2mpZ69zMmXlXJklSUyttgHMSQxV95jNw881w8cXZjg6rrZa1yLn8iCRJuShtgHMSQw1suinceCNcemk2Rm6DDeDWW/OuSpKkplPaAKca2mSTbNbqMcfAgQfCF78Izz2Xd1WSJDUNA5x6JwJ23x2efDJriRs+HMaMgf/+N+/KJEkqvYXyLqAMWlth7Ni8q6iv1tYsrzFgABx9dNYSd/TR2dIjP/gBHHAAtLTkXaYkSaUUKaW8a6iJiBgBjBg2bNgh//jHP/Iup3TGjp1HaH3gATj8cHj77Wyrri23rHNlkiQVV0Q8kFJqm9/zStuF6iSGnGywAdxxBxx1FOy/P+y2m+PjJEmqstIGOOUoAvbYIxsft956sOGGWaCbNi3vyiRJKgUDnGqntTWbqfr3v8NLL2U7Opx3nuvHSZK0gAxwqr3ll8/WjbvmmizADR+eLUMiSZJ6xQCn+mlry8bHjRkD++2XLUPy/PN5VyVJUuGUNsC5lVaDioA994SnnoJPfzoLdUcf7fg4SZIqUNoA5yzUBtfaCscem42Pe/HF98fHzZmTd2WSJDW80gY4FcQKK8CFF8LVV8O552bj4+64I++qJElqaAY4NYbhw+HOO+G734UvfSlbhuSFF/KuSpKkhmSAU+OIgL32ytaPW3vtbFHgY45xfJwkSd0Y4NR4Fl4YjjsOHnkEJk/O9le94ALHx0mS1KG0Ac5ZqCWw4orw29/CH/4AZ5+d7ehw5515VyVJUu5KG+CchVoiG20Ed90F3/oW7LNPtgyJ4+MkSU2stAFOJRORhbennoI118zGxx17LLz1Vt6VSZJUdwY4FcvCC8MJJ2Tj4154IVs/7sILHR8nSWoqBjgV04orwkUXwf/7f/DrX8Pmm2ezVyVJagIGOBXbxhtn4+P22ScLcSeeCDNn5l2VJEk1ZYBT8fXpA1/7Gjz0ENx3H6y/PtxzT95VSZJUMwY4lcdKK8GECdnkhl12gSOOcJKDJKmUDHAql87dHB57DN54A9ZaC264Ie+qJEmqqtIGOBfybXJLLZXNTh0/Hr76VdhvP3j11byrkiSpKkob4FzIVwBsuy08+igMHpy1xl1yCaSUd1WSJC2Q0gY46T0DB8Ipp8A118BPfgI77gj/+lfeVUmS1GsGODWPDTeEiRNh002znRzOPBNmz867KkmSKmaAU3Pp1w+OOQbuuAN+//ts7bgnnsi7KkmSKmKAU3NaYw247bZscsOWW8L3vw/vvJN3VZIk9YgBTs2rT59shupDD8EDD2QLAN99d95VSZI0XwY4acUV4eqr4YQTYNdd4RvfgGnT8q5KkqR5MsBJkC0AvMce8PjjWXhbay24/vq8q5Ikaa4McFJXSy4J558P556b7a+6774wdWreVUmS9AGlDXDuxKAFsvXW2QLAyy0Ha68NF13kAsCSpIZR2gDnTgxaYIssAj//OVx7LfzsZ9kCwC+9lHdVkiSVN8BJVdPWli0A3NYG660HV12Vd0WSpCZngJN6om/fbK24q66C73wHDj7YmaqSpNwY4KRKbLIJPPxw9vl667lunCQpFwY4qVKLLgq/+U02Lm6XXeD442HWrLyrkiQ1EQOc1Fu77JLt4nD//bDZZvDMM3lXJElqEgY4aUEstxxcdx3svz9suimcfbbLjUiSas4AJy2oiGzR3zvugPHjYeRIePnlvKuSJJWYAU6qlk9+MpvUsPbasO66MGFC3hVJkkrKACdVU79+8OMfw+WXwze+AV/5Crz9dt5VSZJKxgAn1cLmm8Mjj8A772TLjdx3X94VSZJKxAAn1cpii8EFF2QtciNGwIknwrvv5l2VJKkEDHBSre22Gzz4INx5Z9YyN2lS3hVJkgrOACfVwworwA03wF57wcYbw7nnutyIJKnXShvgImJERIxvb2/PuxQp06cPHH443Hor/OIXsPvu4O+nJKkXFsq7gFpJKU0AJrS1tR2Sdy3SB6y1FtxzD3z727D++vD730NbW95VLbBx42D69I9+TmsrjBlTn3okqcxKG+CkhjZgAPzyl3DFFbD99tl+qocdli0KXFDTp8PYsR/9nPk9LknqmdJ2oUqFsPvu2eK/F1yQTXZ48828K5IkFYABTsrbsGFw112w/PJZl+r99+ddkSSpwRngpEbQv382seGnP4UddoDTT3eWqiRpngxwUiPZbbesS/V3v4Ndd4U33si7IklSAzLASY3mE5+Av/0NVlop61J1Gy5JUjfOQlUp9GQJiwVV1yUw+veHM86ArbaCHXeEo4/O1pAr8CxVSVL1GOBUCj1ZwmJB5bIExq67wrrrwp57ZgsAn38+LLFEDoVIkhqJXahSo1tllWwf1SFDYL314N57865IkpQzA5xUBP37w2mnwamnwogRcMopzlKVpCZmF6pUJLvskrXC7bkn3HEHXHghLLZYzd6ukrGFra01K0OS1I0BTiqaIUPg9tvhm9+E4cPhqqtgzTVr8lb1GFsoSaqcXahSEfXvD7/6VTY7dcst4fLL865IklRHtsBJRbb//rDOOtls1Xvvzfo8F/KylqSyswVOKrr11oOJE+Hxx2HrreHll/OuSJJUYwY4qQyWWgr+9CfYfHNoa4N77sm7IklSDRngpLJoaYEf/AB++UsYORJ+/WuXGpGkkjLASWUzcmS2l+qvfgUHHlj7PcYkSXVngJPKaNVVs27UmTNhs83g+efzrkiSVEUGOKmsFlkELr44m6m68cZw4415VyRJqpJCBbiIWCUizo2IK/OuRSqECDj8cLjiiqw79Yc/hDlz8q5KkrSAcg9wEXFeRLwSEY91O75dRDwdEZMi4kiAlNJzKaVR+VQqFdgWW2RLjVx3Hey2G0yblndFkqQFkHuAAy4Atut6ICJagF8C2wNrAntHRG32CpKaxfLLw1//CksuCZtsAs8+m3dFkqReyj3ApZRuB17vdnhDYFJHi9tM4DJgp7oXJ5VN//5wzjnw1a/CppvCn/+cd0WSpF7IPcDNwwrA5C73pwArRMRSEXEWsF5EHDWvF0fE6IiYGBETp06dWutapWKJgK99DX7/e9hvPzj1VNeLk6SCadQAF3M5llJKr6WUDk0pfSKldNK8XpxSGp9SaksptQ0ePLiGZUoFttVWcPfdcOGFcMABMGNG3hVJknqoUQPcFGClLvdXBP6dUy1SeQ0Zki36O2NGNtHhxRfzrkiS1AML5V3APNwPrBoRQ4EXgb2AffItSSqpRRaByy6Dn/wENtwQrrwym+SQg3HjPnrjiNZWGDOmfvVIUqPKPcBFxKXAVsDSETEFOCGldG5EHAbcCLQA56WUHq/wvCOAEcOGDat2yVL5RMBRR8E668BOO8FJJ8Go+q/YM306jB0778c/6jFJaia5B7iU0t7zOH4dcN0CnHcCMKGtre2Q3p5Dajo77AC3356FuIcfps/ipwB9865KktRNo46Bk5SXNdaAe++FZ59lv4u2hddey7siSVI3BjhJH7b44jBhAv9erg022giefDLviiRJXZQ2wEXEiIgY397enncpUjG1tHDztj+DY4+FLbeEG2/MuyJJUofSBriU0oSU0uhBgwblXYpUbAccAH/4Q/bxjDNc9FeSGkBpA5ykKvrMZ+Cuu2D8eDj0UJg1K++KJKmpGeAk9czQoVmIe/FF+PznndwgSTkywEnqucUWg6uvhvXXh403hqeeyrsiSWpKpQ1wTmKQaqSlBX7+82zh3y22gJtuyrsiSWo6pQ1wTmKQauygg7Jtt778ZfjFL5zcIEl1lPtODCqm1tbG2taotTXvCprUFltk4+JGjIAnnshmqfZ15wZJqjUDnHrFDcX1nlVWgbvvhr33hu22gyuugCWXzLsqSSo1A5yUk3Hjss3be6u1tYGC9GKLwTXXwHe/m01umDABVl8976okqbQMcFJOpk9fsG7oceN69/qaBb+WFjjlFPjUp2DzzeHii2GbbWrwRpKk0ga4iBgBjBg2bFjepUg10dsQVvOxi6NGwbBhsOeecNxx8LWv1fgNJan5lDbApZQmABPa2toOybsWqelsueX7kxsefxxOPz3XyQ3z6q5uqG5oSapAaQOcpJytskoW4vbeG7bfHi6/HMhncsO8uqsbaSa1JFWitOvASWoAgwZlExrWXhs23pilXn0674okqRQMcJJqq6UFTj0VvvtdDrxgC7jllrwrkqTCM8BJqo9DDuHKL14G++wD48fnXY0kFZoBTlLdvDD0s3DHHdleqt/6FsyenXdJklRIpQ1wbmYvNajVVoN77oGHH4addoJp0/KuSJIKp7QBzs3spQa25JJw442w/PKw2Wbwz3/mXZEkFUppA5ykBte3L5x9Nhx4IGyySdYqJ0nqEdeBk5SfCPjmN2HVVWHECNba4gxg77yrkqSGZwucpPztuCPccgtb33IUnHACpJR3RZLU0AxwkhrDOutwzsH3wk03wV57zX3vK0kSYICT1EDeHrgM/PWv2eK/W20FL72Ud0mS1JAMcJIay4ABcPHFsMMOsNFG2XIjkqQPKO0khogYAYwYNmxY3qVIDaW1teebuLe21rSUeYuA44+H1VeHbbaBc8+FkSNzKkaSGk9pA1xKaQIwoa2t7ZC8a5EayZgxeVdQgT33hCFDYJdd4JlnIH0biLyrkqTc2YUqqbFttFG2RtzvfsfICYfAzJl5VyRJuTPASWp8K68Md97JIm+/AttuC6+9lndFkpQrA5ykYlh0US7b8yoYPhw23hiefjrviiQpNwY4SYWR+rTAz34GRx4JW2wBt9ySd0mSlAsDnKTiGTUKfv972GefbD9VSWoyBjhJxbTVVnDnnXDKKdl+qrNn512RJNVNaZcRkdQEVl01m6G6226w005wySWw2GI9fvnc1sRrbS3YUiuSmpIBTlKxLbEE3HADHHYYbLYZTJiQrR3XA3MLaj1d5FiS8mQXqqTi69sXzjorGxu36aZw9915VyRJNVXaFji30pKaTAQccUTWrTpyJJx2Guy7b95VFdq4cTB9ut3KUiPqUQtcRCwUEXtERP9aF1QtKaUJKaXRgwYNyrsUSfW0ww7wl7/Ascdm+6nOmZN3RYU1fXrWpTx9et6VSOquRwEupfQucG5K6Z0a1yNJC27tteHee+HPf4a99qLvrP/LuyJJC2jcuOymTCVj4O6PiHVqVokkVdPHPpa1xPXrx4Hnbw6TJ+ddkaQFMH26rcFdVTIG7q/AhIgYD/wTeK9fIqV0SbULk6RKdI7X+qAB8InfseUbJ7P8RhvBFVdkM1UlqeAqCXAHkYW2g7sdT4ABTlKuOsdrfVgA34HrPwW77AI//jEc3P2fMUkqlh4HuJTS0FoWIkk1tf32cMcd2QzVhx+GU0/Nlh+RpAKqeBmRiFgeWDmldE8N6pGk2ll99Wxywz77wOc/D5dfDksv3aOXzq2L1uU1VCudg/X9/dK89DjARcTHyLpK/wf4P2BgROwJbJlS+t8a1SdJ1bX44tluDUcfDRtuCFdfnc1anY+5ddG6a4NqxcH6mp9KZqGeATwPDAZmdRz7C7BttYuSpJpqacmaOH7wA/if/4Grrsq7IkmqSCVdqJ8FPp5SmhERCSClNDUiBtemNEmqsX33zbpVd90V/v53OO443GFQUhFU8i/VO3QLfBGxJPB6VSuSpHpqa4P77oMbb4Tdd6ffzLfyrkiS5quSAHcTcHJEdJ22NRb4U1UrkqR6W3ZZ+OtfYfHFGXXupjBpUt4VSeqlZtmxoZIA9z3gk8AbwGIR8SawDnBsLQqTpLrq3x9+8xsmth0Km24K116bd0VS01qQCULNsmNDJevAvQ5sERFtwBCy3RgmppRSjWqTpPqK4P7h/8sOx64He+yRda2ecEI26UFSzTVDy1m19LgFLiI2AEgpTUwpXZlSur+Rw1tEjIiI8e3t7XmXIqloNtkEJk6E22+HHXeE1x3qK9VDs7SeVUMlXah/iYhXI+KKiPhKRHyiZlVVQUppQkpp9KBBg/IuRVIRLbMM3HwzrLkmtLWx7EsP5V2RJL2nkmVElgQ2Bj4H7AucHhEvATenlEbXojhJylXfvnDyybDhhux30Law6c9h//3zrqrhdO5S4c4UUv30uAUupTQ7pfS3lNKJwNeBnwBLAP5rJqnc9tyTC/a/FX78Y/jqV+Gdd/KuqKF07lJh15dUP5WMgds/Ii7qaHW7CFgc2A/o2UaCklRgUz/2qWxSw3/+A1tuCVOm5F2SpCZWyRi484HhwDeAdVNKR3SMM5tWm9IkqcEMGgR/+APsvDMMH87Q527JuyKpoblfcO1UMgZuG7Lxb98FxkfEHcDNZGPgnqpFcZLKpbX1o/9Bb22tWym9FwFHHgnDh7PrzvvB8QfD8cfDQpX8cyqp2pptCZJK1oG7BbgFODoiFgcOA34AnAa4SJKk+SrVAPfPfY6zv/Ig37nrS/C5z8Gll8Lyy+ddldS0mm0MZiVj4JaNiC9FxAXAo8AxwIO4E4OkJvXWwGWzPVS33ho22ABuuCHvkiQ1iUra/KcAj5C1wo0Cbk8pzahJVZJUFC0tcNxxsMUWsO++sN9+cOKJ2RIkklQjlUxiWCaltEFK6XsppZsMb5LUxZZbwkMPZbettoLJk/OuSFKJVTIG7rWIGAjsAKwETAb+lFJ6q1bFSVKhDB4M110HP/0ptLXBb34DI0bkXdUC6zr5xMV6pcbQ4wAXEZ8im3U6G3iBbEP7UyNi25TSYzWpTpKKpk+fbJbq5pvD3nvDX/4CJ50EAwbkXVmvdQ1sLgshNYZKulBPA84GVk4pbQ6sDPwaOL0WhUlSoW22WdadOnkyDB8OjzySd0VSzTVawB83rrzLi1QyiWE94AsppQSQUkoR8RPgiJpUJklFt9RScMUV8NvfZjNVx4yBb30ra6Xr0LmPaCe7KNUIOkNP0X8Xy7y0SCUBrp2s2/QfXY4NAf5bxXokqVwiYP/9s1mq++2XjZG78EJYaSXg/X1EOzVaC4aaU5mDT1lUEuAuBP7U0er2PDAU+B5wQQ3qkqRyGToUbrstm+CwwQZw+unZGLk669riZ2tfcyhLa5o+qJIxcD8i2w91DHAtWXi7sOO4JGl+WlrgqKPg+uuzteL22YcBM96sawmdLX5jx9rK0iymTy/uz7rWY9iK3OLd4wCXUpoN/AIY2+V2Zkrp3VoUJkmltcEG8MADsOSSfPXX67iDg0qnWsGoyOGz1irZSqsNeA44CRgJjAOe7TguSarEwgvDmWdy9chz4atfhQMOgNdfz7sqqS6K2PLVaDVX0oX6K+DklNKQlNLmKaWPAz8nW0pEktQLz31iG3j0UVh0UVh7bdZ46o95l9SwyrwkhFSpSgLcJ4GTux07BVijeuVIUhMaOBB+8Qu49FK2ufl7sNdeMHVq3lXVxLhxWUtGb4KY3WnS+yoJcA8Da3U7tnbH8bqIiEUi4sKIOCci9q3X+0pSXWyxBWcd+nC2xMjaa8Nll0G29GZpdE6iMIhVzhZIdVVJgLsJuDYixkbEgRHxfeAa4KaI2KfzVmkBEXFeRLwSEY91O75dRDwdEZMi4siOw7sCV6aUDiEbhydJpTKr78Lws5/BNdfAD38II0fCCy/kXZYagC2Q6qqSAHcQMAvYHzge+DLwbsfxH3XcftiLGi4Atut6ICJagF8C2wNrAntHxJrAisDkjqfN7sV7SVIxbLhhNlN1442zWas//jG8807eVUlV02iTArrq7OZv5BbPSpYRGdqD2yqVFpBSuh3oPvVqQ2BSSum5lNJM4DJgJ2AKWYj7yNojYnRETIyIiVNLOo5EUhPo3x+OOQYmToS774ZPfxpuuSXvqqRc1StYNXqLZyUtcPW0Au+3tEEW3FYA/gB8MSJ+DUyY14tTSuNTSm0ppbbBgwfXtlJJqrWhQ2HChGwXh1GjYJ994KWX8q6qKTR6K8z8jBuX7bhRJo0erOqlUQNczOVYSim9nVI6MKX01ZTSxXWvSpLyNHIkPP44DBkC66wDp50Gs2blXdUCm9/M1M7H8wgiCxoW8g6A06dXZwutZtqvtyhfX6MGuCnASl3urwj8O6daJKlxLLJINh7u9tuzHRzWWiub8FDg2arzmpnaGdwg+9iTILIgy5TUQtlbi/IOqJ2KErqqqVED3P3AqhExNCL6AXuRzXjtsYgYERHj29vba1KgJOXqk5/MAtzpp2f7q269NTxct1Wd6qIz2FXSgrQgy5Q0ShgpkrIH1EaWe4CLiEuBu4HVI2JKRIzq2F/1MOBG4Eng8pTS45WcN6U0IaU0etCgQdUvWpIaxXbbwSOPwG67ZZ+PGuX4uF4yjKhIcg9wKaW9U0rLpZT6ppRWTCmd23H8upTSaimlT6SUfpR3nZLUsBZaKNtP9emnYemls0WAjzsO3nwz78ok1UjuAU6SVCWDBmV9gBMnwpQpsNpq8JOfwNtv512ZpCpbKO8CJEnva2394IDs1tYPjgEbN+79br7uj71nyBA4/3x46ik44QQYNgyOPBK+8hVgQO2KV6F0jverxixV1V9pA1xEjABGDBs2LO9SJKnHuv8x7T67rnOQ/twe+5A11oDf/z6b3HDccXDyybStfRTMOBAGGOTmpVmCTbOO9yvLz7e0XahOYpCkDuuumy0EfPnlrDrpOlhlFfgGwr2NAAAaRElEQVT5z+n3zrSavF2jLeVRKSczvK+oP8OPUpafb2kDnCSpm4035tK9J8D118PEiRx+xirw/e/D6913M1wwC7KUhypXyzXQ/Bk2rtJ2oUrV1n1sUjXOJ+Xi05+Gyy7jvK8/w9f/NS4bI/elL8E3vpF9rtLr3GLLgFZcBjiph4o+XkLq7rWlVoNfnAsnngi/+hVsuilssgkccQRstRVz39VQZdDZStosOxiMG1e+f8NLG+CcxCCpmXWdrdrVXFt+V1gBfvQjOOYYuOgi+NrXoG9f1l/5MHhrbxg4sOb1Nqt6DqgvY4jpqTK2NJY2wKWUJgAT2traDsm7Fkmqt66zVXts4YVh9Gg45BC4+WZWPfzXsPIY2HPPbAmSddetRalNbX7BopoBr4whphJlC7BOYpAkfVAEbLstv9/zKnj0UVhuORg5EjbaiA0eGA9vvJF3hfPV2lqOGZRlmTHZCLp/H4v++1HaFjhJykP3rst5LrZbFCusAMcfn3WvXn89q4y5EIZ8F7beGvbbD77wBaBf3lV+yJgx7y9nUvifgWqi6MHYACdJVdS967I0g8RbWmDHHbli4o586og34cor4dRT4eCD2WHI7rDtfpA2oZEmPnSGtkp+Bp1jBIv+xz0P82rRKsvCuY2mtF2oETEiIsa3t7fnXYoklcvii8PBB8Ntt8EDD/DfQSvBqFEccfoQOPxwPv7CbTB7dt5V9sqYMQaN3ppX6LUbuDZKG+DciUGS6uDjH+eOzY+GJ57g4n2ug8GD2e7Gb8JyyzHymoPhT3+CGTPyrjIXnWutNYp5jQvs7GrOQy3HoZWm9XseShvgJEl1FMHUj30Kjj2Ws7/yINx3H698bK3sL/Qyy8BOO8GvfsUSbzyXd6V1M316Y7XmjRkz95awPFvHFuS9qx3+ihb4DHCSpOobMoR7Nj4Cbr8dnn0W9toL7r2XUeduCqutxvbXfR2uvRbqMMylLDNSe2rcuMq+3gVtJczr+1vN4Fm08AYGOElSrS29NOy9N1x4ISd/+99w+eVMW2wFOOUUWGEFRo/fAL75TdZ48iqYOrXqbz+vlqeyqnTM2YK2Epb9+9uo4c4AJ0mqis79gj+qRSdFH1h3Xe78zJHwl7/Aa69x/fa/gGWWYYMHz8n2Yl1zTTj0UNZ55Hfw1FMwZ07dvga9r1GDizIuIyJJqopeteT078/klTaFIzfl4hlHMvbYd+Hvf4fbbmO1O6+FL5wAr70G668PbW186oXh8PxwGDIkW3C4QXzUUhmdXYyNNB5OxVfaAOdeqJJUQAstlIW19dfnynZYayxZgJs4Ee6/n7UfvQQ+8014+21Ya60P31g6l7I/qgtxzBhbs1R9pQ1w7oUqSSWx1FLw+c/D5z/PZe92hKGpU+Hxx+Gxx7LbZZfBY4/x7dmt8Le1YI01su7YjlvL7KE04o4RUm+VNsBJkkps8GDYaqvs1iklzvn2i3xrm0fhmWfgH/+A66+HSZM46oXJvHHm8rw5eFWGbj0MPvEJWHllWGml7LbsskBLTl+MFkQzzTDuygAnSXpP171cG2kR2k5dJ0p8aExZBP9dbEXYfkXYfvsPPLTQrFks8c9/8qexkxi65iSYNAn+9jeYMgUmT4bXXuObrcvCzVmg2/ZfK8ISK/Gpx5bl7YHLwBMfy9azW2IJ6OP8v0ZSjxmwjdgFboCTJL2n+16ujaY3+5sC0LcvDBvGpGHD4LC5PD5zJud/90WO2C0LdG+fNxn+8Q/WfPJOBr71MtzzCrz8Mrz1Vtb6t8wyPP/2x3h74DKs9dnBbPHQEkwfsARcumQW8pZYgiVfWxJeXSLbesw/t6oyf6MkSQ3hI1vXqvgec50R2q8fby4xFDYfCsDfnoFtxsIVY7OH3wuM77yTjb97+WX+dsorDHz7ZdZadip9Z73Bov99Ea5+A954A15/nS89+wZc8ga0t3PkQovABUvAootmt4EDYdFF2fnZRZnZbyC88/4xBg5kzcez4zP7DeTdvq3weCtLvD4g+/z1AdkX0r9/XVoD69FF2dv3aOT/bNSaAU6S1BB63bpW4Xss0Pn794cVV4QVV2TSqh2BcA5M/1z28PAu5z5jbMd7zZnDaUf/lyO/8gZMm5a14k2bBtOm8cJFb9HvnWmwWHZs4nUv0/+daaz1f9MYMDt7rM/MGXDndL780gwWenc6nDcjayp95x3o358jGQBnt8KALNiNfrOVdxcaAHcMYN9/9YMH+7LHpL7wVF92ebIfTO7Ljo/0hdf6su3EfjCjL5+9qy+09MtaKvv2ZaN7+8JZ/Vjvwb7MaenL7D594fIW6NOHNZ5sgWtaoKWFYZNa4OYWhj7fQoo+zIkWzh3VwgrRAhNbWPalFlKf7DGezF5DSwuLv9nCnMge4z99iFezzwcQEJE9f1rQb2aQCPi/YKFZHcej49hcPjYTA5wkqe56suhvEfQoEPbpw4wBi8PQxT/00MOPZR+/cFT28dou5+o879ix2e30sR88zpw58M47nDZ2BkcePh1mZMFuwukz6PvudA7aZwb3nj+LVfeYyWMXz2LNnWbx3IyZfHqTWfznP7Ng1Vm89fRMGDSLOS2zYObMLFzOmsWSr82Ch2ex8uSZ9Jkzi5bZs+CK2TBnDus+Phv+bzbMns3Gz8yG/8xmi2dnE2kOfebMJtJs+syZDQ/NZucXs/uR5sBfstcwezYHvD77vedy0RwOm9b5vASk7OMv5vCdmR2fn5o4clbnY3M+8Lwgvfc9SwQpAjo+pujz3ucQ8LPg6FkffIzTgzEzOp5zZvDd6X0+cI7Oj68uvQaMvWX+vxR1UtoA5zpwkhpBZ1Dpet8FXevT2lYrDTPrsU8faG1lRmsrLL/Ee4dfWr7jk61h0p3ATvDEQ8De8MjTsMvBMHEK7PgNuOt12PYouO0d+OzY9099/VjYaCxc3eXYWh2fXzb2/Z/bRR2fX9jleZ3GjoWzxn7wfqfTuh3/6Txe/+Ox73/+w7k85z2dQW5uH9McgkTrgOz+jOkffOzIMYkzfpId+953E7/6aXa8+7nm9FmI73xECfVW2gDnOnCSGkH3sFbEwKIPKvO+n4X1XjcqXdrjPuiY4zquv/7dHlgSpi/c8flgeHuRmlVZVaUNcJIk6cPGjVvwruuGaYVsYi5mI0lSg6tG6Oo0ffqCd+NPn25LZN4McJIkNbhqhC6ViwFOkiRpPhptxrQBTpKkOho3bsHGkC3o69U7jdYC6iQGSeqi+7IfXY9L1bCgY8eKNvascx27RtFo9fSWAU6Suqjkf9lzC3sGPRVdGcJNMzDASVIvNVqXioqr2l2inXu+5skgWFuOgZMkKWfVXpZjzJjidbWqMqUNcBExIiLGt7e3512KJKkCjbpP6rhx79fVeVPt2IL40UrbhepWWpJUTI3aNT19euP/US+ieQU1WxA/Wmlb4CRJKoO8W/uq/d7dQ7BBrXcMcJIk1VklkwzGjKluq2Slgaya720LZvUY4CRJqrM8Jxk0ahe1KlPaMXCS1Ii6rx3nQHhJvWGAk6Q6svVDqpxdrx9mgJMk1dS4ce93FzZai2PnODSDdf7LdqgyBjhJ0nshq5oBq+t6bo3aglKNcWid37Oiz6Ysev3NxgAnSarJGmfN0qrV+XVW+v0rS/Aro0ZrKZ4bZ6FKklRl48bNPwRUe3mQZlLrFt0i/FwMcJIkVdn06cUIASouu1AlqaDmtiSJoUFF1tpa/y7lRh2fOT8GOEkqqO5hrah/iJQp+pi4sWMX/HdwzBh/j3vKLlRJkhqAY+JUidK2wEXECGDEsGHD8i5FklRDlS6B0nV5E6moStsCl1KakFIaPWjQoLxLkSTVUOcSKD3teuzspitqV2VR9WRmrnqutAFOkiTlq2tgc2ZudRngJElSTTRaYCvTBAkDnCSpKbW22qXXU62t7pXaaAxwkqSm5KzPnhszxjGDjaa0s1AlSY3PFjCpdwxwkqTc2ALWGBp5EeE8dmcoArtQJakkOtc3a4Y1zjq/VsdlVUe1u5OrOb6wp3WVaYJCT9gCJ0kl0UytWZ1fa7P90S6Kav8u2gr3YbbASZJUB+PG2WLYW830n5OesgVOkqQ6sAVJ1WQLnCSpqfVkb9Rm2waqmb7WorIFTpLUY13DTlm6tXrydXTut9os8vjZNtP3txoMcJKkHnPygNQY7EKVJEkqGAOcJElqekXb79UAJ0mSml7R9ns1wEmSpIpVc7eFBdGs4zGdxCBJUo018l6jvVWGWcg93eGhEUOiLXCSJNVYtfcaLbt6tewV+WdigJMkqQqKNgi+UY0dW+xgVS+FCnARsUpEnBsRV+ZdiyRJXRVtEHw11aqLsRG7LhtF3QJcRJwXEa9ExGPdjm8XEU9HxKSIOPKjzpFSei6lNKq2lUpSc+jcVWF+20g1skYZSC/VWz0nMVwAnAn8tvNARLQAvwS2AaYA90fENUALcFK31x+UUnqlPqVKUvmVoZuqDF+D1Bt1C3AppdsjYki3wxsCk1JKzwFExGXATimlk4Ad61WbJElSkeQ9Bm4FYHKX+1M6js1VRCwVEWcB60XEUR/xvNERMTEiJk6dOrV61UqSpFIoetd73uvAxVyOpXk9OaX0GnDo/E6aUhoPjAdoa2ub5/kkSb3TOX6u83PpozTi70jRu9/zDnBTgJW63F8R+HdOtUiSeqjof/yaVV6zOqv1+zKvINiMs1XzDnD3A6tGxFDgRWAvYJ98S5Ik6X2N2HrUrPyPw/vqFuAi4lJgK2DpiJgCnJBSOjciDgNuJJt5el5K6fEqvd8IYMSwYcOqcTpJUpMyNKgR1XMW6t7zOH4dcF0N3m8CMKGtre2Qap9bkiQpT3nPQpUkSVKFDHCSJGm+8t7rtRknKnyU0ga4iBgREePb29vzLkWSpMJr5r1eG1FpA1xKaUJKafSgQYPyLkWSpHlyP1f1Rt7LiEiS1NSc5areKG0LnCRJUlkZ4CRJkgqmtF2oLuQrqQy67jnaeV+CfGeEFsnYsR89g7Wo11RpA5wL+UoqA8dHaV6cEVodRb3G7EKVJEkqGAOcJElSwRjgJEmSCsYAJ0mSVDClncTgLFRJUjOpxqzUos7IrKaifA9KG+CchSpJaibVmJVa1BmZ1VSU70FpA5wkSSqvorSU1YoBTpKkghs3Lgs0zbQ23IK0lH3Uwr5F4SQGSZIKbvr04nT9NXvLWbUY4CRJUt0UJWg2OgOcJEkF0drauxasMnQZ6oNKG+AiYkREjG9vb8+7FEmSgN4HsE5jxtiCpUxpJzG4jIgkqdEYvlQtpW2BkyRJKisDnCRJEsUaK1jaLlRJkoqoGZbZ6Pwam2ndumozwEmS1ECaYZxc59dYpBavRmMXqiRJytXcgpzh7qMZ4CRJ6qa1NQsQzdCdqWIqbRdqRIwARgwbNizvUiRJBdMM3Zg9UY9WMFvaeqe0LXAppQkppdGDBg3KuxRJkupiQRcKVnGUtgVOkqRmY8th8yhtC5wkSVJZ2QInqTA6B5bP6zFJahYGOEmFYfeQJGXsQpUkSSoYA5wkSVLBGOAkSZIKxgAnSZJKpRnWwzPASZKkUhkzpvyTnkob4CJiRESMb29vz7sUSZKkqiptgHMrLUmSVFalDXCSJEll5UK+kqRSKPug9SLwZ1A/BjhJUimUfdB6EfgzqB+7UCVJkgrGACdJklQwBjhJkqSCMcBJkiQVjAFOkqQqaW2FsWPLOxuzGbaoKgpnoUqSVCVln4XZ+fWNHVud81XrPM3IFjhJkgrK1rDmZYCTJKmgyt7it6AWtIWvkVsIDXCSJEkFY4CTJEkqmNIGuIgYERHj29vb8y5FkiSpqkob4FJKE1JKowcNGpR3KZIkSVVV2gAnSZJUVgY4SZKkgjHASZIkFYwBTpIkVaSR10drFm6lJUlSCblLQ7kZ4CRJKiF3aXhfGVsM7UKVJEkqGFvgJElSKTRTt7EBTpIklUIzdRvbhSpJkjQPjTp+zgAnSZJUMAY4SZKkgjHASZIkFYwBTpIkqWAMcJIkSQVjgJMkSSoYA5wkSVLBuJCvJEkF00w7DmjuDHCSJBVMM+04oLmzC1WSJKlgChXgImLniDgnIq6OiG3zrkeSJCkPdQtwEXFeRLwSEY91O75dRDwdEZMi4siPOkdK6Y8ppUOAA4A9a1iuJElSw6rnGLgLgDOB33YeiIgW4JfANsAU4P6IuAZoAU7q9vqDUkqvdHx+bMfrJEmSmk7dAlxK6faIGNLt8IbApJTScwARcRmwU0rpJGDH7ueIiAB+AlyfUnpwXu8VEaOB0QArr7xyVeqXJElqFHmPgVsBmNzl/pSOY/PydWBrYLeIOHReT0opjU8ptaWU2gYPHlydSiVJkhpE3suIxFyOpXk9OaV0BnBG7cqRJElqfHm3wE0BVupyf0Xg3znVIkmSVAh5B7j7gVUjYmhE9AP2Aq6pxokjYkREjG9vb6/G6SRJKp2xY/OuQL1Vz2VELgXuBlaPiCkRMSql9C5wGHAj8CRweUrp8Wq8X0ppQkpp9KBBg6pxOkmSGlJrq1trNaNIaZ5Dzkqhra0tTZw4Me8yJEmS5isiHkgptc3veXl3oUqSJKlCBjhJkqSCKW2AcxKDJEkqq9IGOCcxSJKksiptgJMkSSorA5wkSVLBGOAkSZIKprQBzkkMkiSprEob4JzEIEmSyqq0AU6SJKmsDHCSJEkFY4CTJEkqmNIGOCcxSJKksiptgHMSgyRJKqvSBjhJkqSyMsBJkiQVjAFOkiSpYAxwkiRJBRMppbxrqKmImAr8s4qnHATUY2prtd9nQc/Xm9dX+pqePr+nz1saeLWC9y+Dev1+9kQRr5VqnMtrpRi8VvI9V29fn+e1Uq/r5OMppcHzfVZKyVsFN2B8Ed9nQc/Xm9dX+pqePr+C502s9+9H3rd6/X42Ui3VfJ9qnMtrpRg3r5V8z9Xb1+d5rTTadWIXauUmFPR9FvR8vXl9pa/p6fPr9TMookb63hTxWqnGubxWiqGRvjfNeK309vVeKx1K34Wq5hURE1NKbXnXITU6rxVp/hrtOrEFTmU2Pu8CpILwWpHmr6GuE1vgJEmSCsYWOEmSpIIxwEmSJBWMAU6SJKlgDHBqShGxSEQ8EBE75l2L1Kgi4pMRcVZEXBkRX827HqlRRcTOEXFORFwdEdvW4z0NcCqUiDgvIl6JiMe6Hd8uIp6OiEkRcWQPTjUGuLw2VUr5q8a1klJ6MqV0KLAH0DDLJ0jVVKVr5Y8ppUOAA4A9a1ju+/U5C1VFEhFbAG8Bv00prdVxrAV4BtgGmALcD+wNtAAndTvFQcA6ZFuiDABeTSldW5/qpfqpxrWSUnolIkYCRwJnppQuqVf9Ur1U61rpeN3JwMUppQdrXfdCtX4DqZpSSrdHxJBuhzcEJqWUngOIiMuAnVJKJwEf6iKNiM8CiwBrAtMj4rqU0pyaFi7VWTWulY7zXANcExF/AgxwKp0q/V0J4CfA9fUIb2CAUzmsAEzucn8KsNG8npxSOgYgIg4ga4EzvKlZVHStRMRWwK5Af+C6mlYmNZaKrhXg68DWwKCIGJZSOquWxYEBTuUQczk237EBKaULql+K1NAqulZSSrcCt9aqGKmBVXqtnAGcUbtyPsxJDCqDKcBKXe6vCPw7p1qkRua1IvVMw18rBjiVwf3AqhExNCL6AXsB1+Rck9SIvFaknmn4a8UAp0KJiEuBu4HVI2JKRIxKKb0LHAbcCDwJXJ5SejzPOqW8ea1IPVPUa8VlRCRJkgrGFjhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBWOAkyRJKhgDnCRJUsEY4CQVUkQsGRE3RkR7RDyQdz2SVE/uhSqpqA4FBgJLdSy6KUlNwxY4SUW1CvDk3MJbRPTNoR5JqhsDnKTCiYgJwP7A/hHxVkTcFhHvRsR+EfEc8HrH8xaOiJ9HxPMR8XpE3BARw7qcZ9GIuLDjsX9GxP4d59mq4/GxEfHnbu99a0Qc2+X+Wh1dua9GxL8i4qTOABkRQyIiddT1RERMi4ibImK5Lq8f2FHjcx2PPx4Rn4mI7SNiasc+jF3rfSsiNq/Nd1ZSURjgJBVOSmkEcDFwYUppIHAC0AJsD6wHLNPx1N8AawAbA8sC9wLXdmmhOw1YFVgTWAfYqeM8PRIRHwNuA/4ALA9sAmwDHNXtqXsCWwArAIsAJ3Z57FxgI+BzwGLAzsB/yPZgfLujpk57A5NTSnf0tEZJ5WSAk1QmR6aU2lNK/xcRS5MFnv9NKb2cUpoJfB9YDtgoIvoA+wLHpZT+k1JqB8ZU+H5fBh5JKZ2dUpqZUnoROKnjeFffTym9mlL6L3AJ0AbvBcA9gENTSs+nzD9SSpNSSnPIAuioLucZ1XFMUpNzEoOkspgDTO5yf2jHx79HRNfn9QVWAgYD/YEXujz2fIXvORTYLCLe7HIs+HAr3ktdPn8bWLTj8yEdH5+Zx/nPBY6LiJXJWufWBXaosEZJJWSAk1QWKaWUutz/Z8fHVVNKU7s/uaMFbiZZiHq24/DQbk97i6zLs6vlu73Hn1NKvQ1VL3TWCDzR/cGU0ksR8SfgQGAJ4I8ppVd7+V6SSsQuVEmllFJ6hay78lcRsQJARCweEbtExMCOLspLgO9HxDIRsRhZ92dXE4H1I2KDiFgoIg7jgyHvt0BbRBwUEQMiok9ErBIR21VQ45UdNQ6JzLCuEy2A8cBBwJeAcyr/TkgqIwOcpDI7BHgauDUipgGPArsDnS11h5N1mz7V8dgEYHbni1NKtwInAzeQdYMuA/yty+P/AT5LNvHgBeAN4CqyJU566iDgYbLJENOAq8kmXHS6iax7uB24pYLzSiqx+GCPgyQ1t4h4F9i6I7w1hIi4FbgppfTjvGuR1BgcAydJDSwitgCGk7UcShJggJOkhhUR9wPDgK/PbSKGpOZlF6okSVLBOIlBkiSpYAxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGACdJklQw/x8NS9QHesYk5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "perdata00 = pd.read_csv(\"perlist00.csv\")\n",
    "f = perdata00['f']\n",
    "per = perdata00['per']\n",
    "\n",
    "alpha_L = 1.0\n",
    "A,f_b,alpha_H,poisson = m.values[0],m.values[1],m.values[2],m.values[3]\n",
    "\n",
    "model = []\n",
    "f_length = len(f)\n",
    "for i in range(f_length):\n",
    "    model.append(((f[i]**(-alpha_L))/(1+(f[i]/f_b)**(alpha_H-alpha_L)))*A+poisson)\n",
    "    \n",
    "plt.figure(figsize=(10,8))\n",
    "plt.loglog()\n",
    "plt.step(f, per, color=\"b\", alpha=0.5, linewidth=1)\n",
    "plt.plot(f, model, color=\"r\", linewidth=1)\n",
    "plt.xlabel(\"frequency\",fontsize=13)\n",
    "plt.ylabel(\"power\",fontsize=13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
